data <- readRDS('data_sim_p_1000.rds')
for(x in names(data)) assign(x, data[[x]], envir = globalenv())
U_full <- data$U
rm(data)

S.data.time <- S.data$obs
S.data.time.log.sum <- sum(log(S.data.time))
model <- SAEM_model( 
  function(sigma2, ...) -n/(2*sigma2),
  function(phi1, phi2, phi3, ...) mean((Y - get_obs(phi1, phi2, phi3) )^2 ), 'sigma2',
  
  # === Variable Latente === #
  latent_vars = list(
    # === Non linear model === #
    latent_variable('phi', dim = G, size = 3, prior = list(mean = 'mu', variance = 'omega2'),
                    add_on = c('zeta(phi1 = phi1, phi2 = phi2, phi3 = phi3, ...)' )),
    
    # === S.data model === #
    latent_variable('b', prior = list(mean = 'barb', variance.hyper = 'sigma2_b'),
                    add_on = c('zeta(b = b, ...) +',
                               'sum(h$eval(b = b, ..., i = c(1,2)))' )),
    latent_variable('alpha', prior = list(mean = 'baralpha', variance.hyper = 'sigma2_alpha'),
                    add_on = c('zeta(alpha = alpha, ...) +',
                               'alpha*h$eval(alpha = alpha,..., i = 3)'))
  ),

  # === Paramètre de regression === #
  regression.parameter = list(
    regression_parameter('beta', 1, function(...) SPGD(5, theta0 = beta,
                                                      step = 0.05, lambda = 1/sqrt(N),
                                                      normalized.grad = T,
                                                      zeta.der.B, N, zeta.B, 
                                                      Z$alpha,  Z$phi1, Z$phi2, Z$phi3,Z$b) )
  )
)
# ---  Initialisation des paramètres --- #
parameter0 <- parameter %>% sapply(function(x) x* runif(1, 1.1,1.4))

parameter0$barb <- 10
parameter0$sigma2 <- 1
parameter0$baralpha <- 1
#===============================================#
load.SAEM(model)
oracle <- maximisation(1, do.call(S$eval, var.true), parameter, var.true)
#==============================================================================#

init.options <- list(x0 = list(phi = c(1,80,4), b = parameter0$barb, alpha = parameter0$baralpha), 
                     sd = list(phi = c(.05, 1.5, .5), b = 1, alpha = .1) )

SAEM.options <- list(niter = 200, sim.iter = 5, burnin = 190, 
                adptative.sd = 0.6)

saem

p <- 4
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 45min 44sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 45.0000
Estimated value 0.0025 0.8979 90.0274 4.9032 0.0040 34.2305 0.5823 18.2806 30.2388
Rrmse 0.0101 0.0059 0.0008 0.0209 0.0098 0.0470 0.1619 0.3906 0.3280
p <- 20
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 46min 29sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 45.0000
Estimated value 0.0025 0.9004 90.1263 4.9105 0.0044 34.6121 0.5214 17.7837 30.0382
Rrmse 0.0136 0.0031 0.0019 0.0194 0.1135 0.0364 0.2495 0.4072 0.3325
p <- 100
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 46min 12sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 45.0000
Estimated value 0.0026 0.9025 90.1870 4.9241 0.0047 34.7273 0.5773 15.3291 27.0420
Rrmse 0.0189 0.0007 0.0026 0.0167 0.1872 0.0331 0.1690 0.4890 0.3991
p <- 1000
U <- U_full[,1:p] #Mise à jours des cov
parameter$beta <- parameter$beta[1:p] #Mise à jours de beta
parameter0$beta <- runif(p, min = -1, max = 1) #initialisation de beta

SAEM.options$RDS <- paste0(params$rds_filename, '_', p)
res <- run(model, parameter0, init.options, SAEM.options, verbatim = 3)
## [1] "SAEM execution time = 36min 26sec"
## $plot_parameter

## 
## $plot_MCMC

## 
## $plot_acceptation

## [[1]]

## 
## [[2]]

## 
## [[3]]
Result of the SAEM-MCMC
sigma2 mu.1 mu.2 mu.3 omega2.1 omega2.2 omega2.3 barb baralpha
Real value 0.0025 0.9032 89.9575 5.0079 0.0039 35.9179 0.6948 30.0000 45.0000
Estimated value 0.0025 0.9037 90.0198 5.0017 0.0040 35.3809 0.6134 4.2174 8.2451
Rrmse 0.0046 0.0005 0.0007 0.0012 0.0277 0.0150 0.1172 0.8594 0.8168